Optimal waveform estimation for classical and quantum systems via time-symmetric 

smoothing 
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Classical and quantum theories of time-symmetric smoothing, which can be used to optimally 
estimate waveforms in classical and quantum systems, are derived using a discrete-time approach, 
and the similarities between the two theories are emphasized. Application of the quantum theory to 
homodyne phase-locked loop design for phase estimation with narrowband squeezed optical beams 
is studied. The relation between the proposed theory and Aharonov et al.'s weak value theory is 
also explored. 
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FIG. 1: (Color online). Four classes of estimation problems, 
depending on the observation time interval relative to r, the 
time at which the signal is to be estimated. 



Estimation theory is the science of determining the 
state of a system, such as a dice, an aircraft, or the 
weather in Boston, from noisy observations P, [3j S 
As shown in Fig. [2 estimation problems can be classified 
into four classes, namely, prediction, filtering, retrodic- 
tion, and smoothing. For applications that do not re- 
quire real-time data, such as sensing and communication, 
smoothing is the most accurate estimation technique. 

I have recently proposed a time-symmetric quantum 
theory of smoothing, which allows one to optimally esti- 
mate classical diffusive Markov random processes, such 
as gravitational waves or magnetic fields, coupled to a 
quantum system, such as a quantum mechanical oscilla- 
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tor or an atomic spin ensemble, under continuous mea- 
surements Q. In this paper, I shall demonstrate in 
more detail the derivation of this theory using a discrete- 
time approach, and how it closely parallels the classical 
time-symmetric smoothing theory proposed by Pardoux 
0. I shall apply the theory to the design of homodyne 
phase-locked loops (PLL) for narrowband squeezed opti- 
cal beams, as previously considered by Berry and Wise- 
man I shall show that their approach can be re- 
garded as a special case of my theory, and discuss how 
their results can be generalized and improved. I shall also 
discuss the weak value theory proposed by Aharonov et 
al. ^ in relation with the smoothing theory, and how 
their theory may be regarded as a smoothing theory for 
quantum degrees of freedom. In particular, the smooth- 
ing quasiprobability distribution proposed in Ref. h\ is 
shown to naturally arise from the statistics of weak po- 
sition and momentum measurements. 

This paper is organized as follows: In Sec. [Ill Par- 
doux's classical time-symmetric smoothing theory is de- 
rived using a discrete-time approach, which is then gen- 
eralized to the quantum regime for hybrid classical- 
quantum smoothing in Sec. IIIII Application of the hy- 
brid classical-quantum smoothing theory to PLL design 
is studied in Sec. IIVI The relation between the smooth- 
ing theory and Aharonov et aZ.'s weak value theory is 
then discussed in Sec. |Vl Sec. IVII concludes the paper 
and points out some possible extensions of the proposed 
theory. 



II. CLASSICAL SMOOTHING 

A. Problem statement 

Consider the classical smoothing problem depicted in 
Fig.m Let 

Xit 
X2t 



Xt 



(2.1) 
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FIG. 2: (Color online). The classical smoothing problem. 



be a vectoral diffusive Markov random process that sat- 
isfies the system Ito diff'erential equation [l| 



dxt = A{xt,t)dt + B{xt,t)dWt, 



(2.2) 



where dWt is a vectoral Wiener increment with mean and 
covariance matrix given by 



(dWt) = 0, 
{dWtdWf^} = Q{t)dt. 



(2.3) 
(2.4) 



The superscript ^ denotes the transpose. The vectoral 
observation process dyt satisfies the observation Ito equa- 
tion 



dyt = C{xt,t)dt + dVt, 



(2.5) 



where dVt is another vectoral Wiener increment with 
mean and covariance matrix given by 



{dVt) = 0, 
(dVtdVt^) = R{t)dt. 



(2.6) 
(2.7) 



For generality and later purpose, dWt and dVt are as- 
sumed to be correlated, with covariance 

(dWtdVt^) = S{t)dt. (2.8) 

Define the observation record in the time interval [ti, 



dy[uj2) = {dyt,ti <t<t2}. 



(2.9) 



The goal of smoothing is to calculate the conditional 
probability density of Xt, given the observation record 
dy[to,T) in the time interval to < t < T . 

It is more intuitive to consider the problem in discrete 
time first. The discrete-time system and observation 
equations (|2.2p and (|2.5p are 

dxt ^ A{xt,t)St + B{xt,t)dWt, (2.10) 
Syt^C{xt,t)St + SVt. (2.11) 

The observation record 

Sy[to.T^st] = {SytoJyto+st, • ■ • , Syrst} (2.12) 

also becomes discrete. The covariance matrices for the 
increments are 



(dWtSW^) ^ Q{t)dt, 
(dVtSVt^) = R{t)dt, 
(SWtSVt^) = Sit)6t, 



(2.13) 
(2.14) 
(2.15) 



and the increments at different times are independent of 
one another. Because 6Wt and 6Vt are proportional to 
one should keep all linear and quadratic terms of 
the Wiener increments in an equation according to Ito 
calculus when taking the continuous time limit. 

With correlated SWt and 6Vt, it is preferable, for tech- 
nical reasons, to rewrite the system equation (|2.10|) as 
i 



Sxt = A{xt,t)St + B{xt,t)6Wt 

+ D{xt,t) [5yt~C{xt,t)5t-5Vt] 



(2.16) 



where D{xt,t) can be arbitrarily set because the expres- 
sion in square brackets is zero. The system equation be- 
comes 

Sxt = [A{xt,t) - Dixt,t)C{xt,t)] St + D{xt,t)Syt 

+ B{xt,t)SWt-D{xt,t)SVt. (2.17) 
The new system noise is 

6Zt = Bixt,t)5Wt - D{xt,t)5Vt, (2.18) 
(SZtSZ^) = [B{xt,t)Q{t)B^{xt,t) 
+ D{xt,t)R{t)D^{xt,t) 
-B{xt,t)Sit)D^ixt,t) 
- D{xt,t)S^it)B^{xt,t)]5t. (2.19) 

The covariance between the new system noise SZt and 
the observation noise 6Vt is 

(SZtSVt^) ^ [B{xt,t)S{t) - D{xt,t)R{t)] St, (2.20) 

and can be made to vanish if one lets 

D{xt,t) = B{xt,t)S{t)R-\t). (2.21) 

The new equivalent system and observation model is then 

5xt = A{xt,t)5t -I- B{xt,t)S{t)R-\t) [5yt - C{xt,t)5t] 

+ B{xt,t)6Ut, (2.22) 

Syt^Cixt,t)St + dVt, (2.23) 

with covariances 

(SUtSU^) = [Q{t) - S{t)R-^S^it)] 6t, (2.24) 

(SVtdVt^) = R{t)St, (2.25) 

(SUtSVt^) = 0. (2.26) 

The new system and observation noises are now indepen- 
dent, but note that Sxt becomes dependent on Syt- 

B. Time-symmetric approach 

According to the Bayes theorem, the smoothing prob- 
ability density for Xt can be expressed as 

P{Xr\Sy[to^T-St]) = ^77 ■ T , 2.27) 

P(Sy[to.T-st]) 

P{5y[to,T^St]) =y dXrP{5y[ta,T-St]\xT)P{Xr), 

(2.28) 
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where 



dxr = / dxir 



(2.29) 



and P{xr) is the a priori probabihty density, which rep- 
resents one's knowledge of Xr absent any observation. 
Functions of assumed to also depend implicitly on 

T. Splitting 5y[to.T-5t] into the past record 



^J/past = ^y[to.T-5t] 

and the future record 

(^J/futuro = ^y[T,T-St] 



(2.31) 



relative to time r, P{5y\^ta,T)\xT) in Eq. (|2.27p can be 
rewritten as 

P{^y[to,T~Si\\xT) = ^'('52/past,%uturo|a;r) 

= P{5yfntuTc\5yvi>i^t,Xr)P{5yp^t\xT)- 

(2.32) 

Because 5Vt are independent increments, the future 
record is independent of the past record given Xt, and 

P{^y[to,T-St\\xT) = P{hintnvc\xT)P{5yp?^st\xT)- (2.33) 

Equation (I2.27|) becomes 

P{hhxtnTc\xT)P{5yp!,st\xT)P{Xr) 



P{Xr\5y[t„,T-5t]) 



J dxr (numerator) 

P{SyfutuYc\xT)P{xT\Sypi>st) 
J dxr (numerator) 



(2.34) 



Thus, the smoothing density can be obtained by combin- 
ing the filtering probability density P{xr\Sypast) and a 
retrodictive likelihood function P{Sy{uture\xT)- 



C. Filtering 

To derive an equation for the filtering probability den- 
sity P{xr\Sypast), first express P{xt+st\Sy[to,t]) in terms 
of P{xt\Sy[t„,t]) as 



Pixt+st\Syito,t]) = / dxtP{xt+st,Xt\6y[to.t]) 



= J dxtPixt+st\xt,Sy[t^^t])P{xt\5y[to,t])- 

(2.35) 

P{xt+st\xt,Sy[to,t]) = P{xt+st\xt,5yt, ^y[to t~st]) can be 
determined from the system equation (|2.22p and is equal 
to P{xt+st\xt,Syt), due to the Markovian nature of the 
system process. So 



P{xt+st\Sy[ta,t]) ^ I dxtPixt+5t\xt,Syt)Pixt\Sy[to,t]), 

(2.36) 



which is a generalized Chapman-Kolmogorov equation 
P{xt+st\xt,Syt) is 

P{xt+5t\xt,Syt) oc 

expj - [B{xt,t)Q{t)B^{xt,t)6tyUZty 



(2.37) 



(2.30) where 



SZt = xt+st ~ Xt ~ A{xt,t)St 

+ B{xt,t)S{t)R-\t) [Syt - C{xt,t)St] . (2.38) 

Next, write P{xt\Sy[tg^t]) in terms of P{xf \6y[to^t-St]) us- 
ing the Bayes theorem as 

P{xt\Sy[t„^t]) = P{xt\Sy[t„,t-5t],Syt) 

P{6yt\xt,5y[to,t-st])P{xt\Sy[to,t-st]) 



J dxt (numerator) 

P{6yt\xt)P{xt\Sy[to^t~st]) 
J dxt (numerator) ' 



(2.39) 



where P{Syt\xt, Sy[to,t-5t]) = P{Syt\xt) due to the 
Markovian property of the observation process. 
P{Syt\xt) is determined by the observation equation 
(|2.23p and given by 

P{Syt\xt) « exp I - i [Syt - C{xt,t)6tf [R{t)Str' 

x[dyt-C{xt,t)6t]}. (2.40) 



Hence, starting with the a priori probability density 
P{xtg), one can solve for P{xr\Sypast) by iterating the 
formula 



(2.41) 



P{xt+st\Sy[to.t]) ^ I dxtP{xt+st\xt,Syt) 

P{5yt\xt)P{xt\Sy[to,t-St]) 
J dxt (numerator) 



To obtain a stochastic differential equation for the filter- 
ing probability density, defined as 



F{x,t) = P{xt = x\dy\to,t)) 



(2.42) 



in the continuous time limit, one should expand 
Eq. (|2.4ip to first order with respect to 5t and second 
order with respect to 5yt in a Taylor series, then ap- 
ply the rules of Ito calculus. The result is the Kushner- 
Stratonovich (KS) equation [H, [l3|, generalized for cor- 
related system and observation noises by Fujisaki et al. 



[ll|, given by 



where 



dt X ^ 
T 2^ 



dx„dx„ 



1 F 

/ uv 



+ {C^ {C)Ff R-^df]tF 



where 



dF = F{x,t + dt)''F{x,t), 
{C)f= J dxC{x,t)F{x,t), 
drjt = dyt - dt{C')F- 
The initial condition is 



(2.43) 

(2.44) 
(2.45) 
(2.46) 

(2.47) 



drjt is called the innovation process and is also a Wiener 
increment with covariance matrix R{t)dt (ill, [l^. 

A linear stochastic equation for an unnormalized F 
is called the Duncan-Mortensen-Zakai (DMZ) equation 

Sin, 

given by 



where the normalization is 



F(x,t) 



f{x,t) 
J dxf{x,t)' 



(2.48) 



(2.49) 



D. Retrodiction and smoothing 

To solve for the retrodictive likelihood function 
-P('5j/futuro|a;r), note that 

^'('^yfuturc) = J dXrPiSy{ntuvc\Xr)P{XT): (2.50) 

but P(5yfuturo) can also be expressed in terms of the mul- 
titime probability density as 



X[r,T] = {xT,Xr+St,---,XT}, (2.52) 

J Dx[t t]= J dxr J dxr+st ■ ■ ■ J dxT ■ (2.53) 
The multitime density can be rewritten as 

Pix[r,T],Sy[r,T-St]) = P{xT\x[r,T-St], Sy[r^T-5t]) 

-X Pixir,T-st],Syir.T-st])- (2-54) 

Again using the Markovian property of the system pro- 
cess, 

P{xT\x[r.T-St], Sy[T.T-5t]) = P{xT\xT-5t, SyT-St) , 

(2.55) 

which can be determined from the system equation 
(1^?^ and is given by Eq. (|07|) . Furthermore, 
P{x[r.T~st]T Sy[r^T-St]) Eq. (|2.54p can be expressed as 

Pixir,T~St],SyiT.T-St]) = P{SyT^St\xir.T~St],Syir,T~2St]) 

X P{xir.T-5t],SyiT,T-2St])- 

(2.56) 

Using the Markovian property of the observation process, 

P{SyT-St\x[r,T~5t],^y[T.T^25t]) = P{SyT-St\xT-St) , 

(2.57) 

which can be determined from the observation equation 
(12:231) and is given by E q. ((2:401) . Applying Eqs. (|234l) . 
(|235l) . (12361) . and ([2371) repeatedly, one obtains 

PiSy[T,T-St]) ^ j dxT J dxT-6tPixT\xT-5t,5yT-5t) 

X P{SyT~st\xT~st) 

X J dxT~2StP{xT-5t\xT-25t,SyT-2St) 
X P{SyT-2St\xT-2St) ■ • ■ 
X j dXrP{Xr+St\Xr,Syr) 

X P{Syr\Xr)P{Xr). (2.58) 



P{Sy[T.T-5t]) = / P'X[r,T]P{x[T.T],Sy[T.,T-St]), (2.51) 



Comparing this equation with Eq. (|2.50p . ^((Jyfuturola^T) 
can be expressed as 

Pi^yiuturelXr) = P{Syr\Xr) 

X J dxr+stP{xT+st\xr,5yr) ■ ■ . 

X P{SyT-25t\xT-2St) 

X J dxT-StP{xT-St\xT-25t,SyT-2St) 
X P{SyT-5t\xT-5t) 

X J dxTP{xT\xT-5t,SyT-5t)- (2.59) 
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Defining the unnormalized retrodictive likelihood func- 
tion at time t as 



gix,t) oc P{dy[t^T)\xt = x), 



(2.60) 



one can derive a linear backward stochastic differential 
equation for g{x,t) by applying Ito calculus backward in 
time to Eq. (f239l) . The result is d 



C^R-Hytg + Y,{BSR-^dyt)^^. (2.61 



dx^ 



must be the same for any r, and one can simply esti- 
mate x^r at the end of the observation interval (r = T) 
using filtering alone. This also means that smoothing is 
not needed when one only needs to detect the presence 
of a signal in detection problems since the presence 
can be regarded as a constant binary parameter within 
a certain time interval. In general, however, smoothing 
can be significantly more accurate than filtering for the 
estimation of a fluctuating random process in the middle 
of the observation interval. Another reason for modeling 
unknown signals as random processes is robustness, as 
introducing fictitious system noise can improve the esti- 
mation accuracy when there are modeling errors P, 0] ■ 



which is the adjoint equation of the forward DMZ equa- 
tion (|2.48p . to be solved backward in time in the back- 
ward Ito sense, defined by 

-dg = g{x,t^dt)- g{x,t), (2.62) 

with the final condition 

g{x,T)ozl. (2.63) 

The adjoint equation with respect to a linear differential 
equation 



is defined as 



df{x,t)^Lf{x,t) 



-dg{x,t) = Dg{x,t), 



(2.64) 



(2.65) 



where L is a linear operator and U is the adjoint of L, 
defined by 

{g{x),Lf{x))^{L^g{x)J{x)) (2.66) 
with respect to the inner product 

{g{x)J[x)) = f dxg{x)f{x). (2.67) 



After solving Eq. ^(IWji for /(a;,r) and Eq. (P3T|l for 
g{x,T), the smoothing probability density is 



h{x,T) = P{xr = x\dy\^to,T)) 



g{x,T)f{x,T) 
J dxg{x,T)f{x,T) ' 

(2.68) 



Since f{x,T) and g{x,T) are solutions of adjoint equa- 
tions, their inner product, which appears as the denomi- 
nator of Eq. (|2.68p . is constant in time The denomi- 
nator also ensures that h{x, r) is normalized, and /(x, t) 
and g{x,T) need not be normalized separately. 

The estimation errors depend crucially on the statistics 
oixt- If any component of Xt, say x^t, is constant in time, 
then filtering of that particular component is as accurate 
as smoothing, for the simple reason that P{Xfj_T\dy[t„,T)) 



E. Linear time-symmetric smoothing 

If /, g, and h are Gaussian, one can just solve for 
their means and covariance matrices, which completely 
determine the probability densities. This is the case when 
the a priori probability density P(xtg) is Gaussian, and 



A{xt,t) = J{t)xt, 
B{xt,t) = B{t), 
C{xt,t) ^ K{t)xt. 



(2.69) 
(2.70) 
(2.71) 



The means and covariance matrices of /, g, and h 
can then be solved using the linear Mayne-Fraser-Potter 
(MFP) smoother [ij. The smoother first solves for the 
mean x' and covariance matrix S of / using the Kalman 
filter [l| , given by 

dx' = Jx'dt + T{dy- Kx'dt) , (2.72) 
T = {Y.K'^ + BS) R~^, (2.73) 

dE = ( js + E - ri?r^ + bqb"^) dt, (2.74) 

with the initial conditions at to determined from P{xtg). 
The mean x" and covariance matrix S of g are then solved 
using a backward Kalman filter, 

-dx" = -Jx"dt + T{dy - Kx"dt), (2.75) 
r={EK^ + BS)R-^, (2.76) 
-dE = (- JS - EJ'^ - TRT^ + BQB'^) dt, (2.77) 

with the final condition S^^x^ = and — 0. In 
practice, the information filter formalism should be used 
to solve the backward filter, in order to avoid dealing 
with the infinite covariance matrix at T [mi. Finally, 
the smoothing mean Xr and covariance matrix 11,- are 

Xr = Ilr (£-14 + S-14') , (2.78) 

n.= (E;i+S;i)-\ (2.79) 



Note that 



and E are the mean and covariance ma- 



trix of a likelihood function P{dy[t.T)\xt) and not those 
of a conditional probability density P{xt\dy\j^T))i so to 
perform optimal retrodiction (r = to) one should still 
combine x" and E with the a priori values . 
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III. 



HYBRID CLASSICAL-QUANTUM 
SMOOTHING 

A. Problem statement 



Smoothing Estimation 



Measurements 

M(dyt\xt) ' 



Classical System 



Quantum System 



dyt 



Observations 



FIG. 3: (Color online). Schematic of the hybrid classical- 
quantum smoothing problem. 

Consider the problem of waveform estimation in a hy- 
brid classical-quantum system depicted in Fig. O The 
classical system produces a vectoral classical diffusive 
Markov random process Xt, which obeys Eq. ()2.2p and 
is coupled to the quantum system. The goal is to esti- 
mate Xr via continuous measurements of both systems. 
This setup is slightly more general than that considered 
in 0; here the observations can also depend on xt- This 
allows one to apply the theory to PLL design for squeezed 
beams, as considered by Berry and Wiseman ^Ti], and po- 
tentially to other quantum estimation problems as well 
[l^ . The statistics of Xt are assumed to be unperturbed 
by the coupling to the quantum system, in order to avoid 
the nontrivial issue of quantum backaction on classical 
systems [13] ■ For simplicity, in this section we neglect the 
possibility that the system noise driving the classical sys- 
tem is correlated with the observation noise, although the 
noise driving the quantum system can still be correlated 
with the observation noise due to quantum measurement 
backaction. Just as in the classical smoothing problem, 
the hybrid smoothing problem is solved by calculating 
the smoothing probability density P{xT\dy[tg T))- 



B. Time-symmetric approach 



respectively. The hybrid operator can also be regarded 
as a special case of the quantum density operator, when 
certain degrees of freedom are approximated as classi- 
cal. Unfortunately, the density operator in conventional 
predictive quantum theory can only be conditioned upon 
past observations and not future ones, so it cannot be 
used as a quantum version of the smoothing probability 
density. 

The classical time-symmetric smoothing theory, as a 
combination of prediction and retrodiction, offers an im- 
portant clue to how one can circumvent the difficulty of 
defining the smoothing quantum state. Again casting 
the problem in discrete time, and defining a hybrid effect 
operator as E{6y{utiiTe\xT), which can be used to deter- 
mine the statistics of future observations given a density 
operator at r, 

P {Sy future) = J dXr tr E {Syiutui-c\Xr) pi^r) , (3.3) 

one may write, in analogy with Eq. (|2.34p 

|. ^ -P(a;T, %uturc|^ypast) 

PiXr\dy[to,T-St]) = |T ^ 

-P(£'2/futuro|02/past) 

_ tr[E{Sy{uturc\xT)p{xT\Sypast)] 

J dXr tr[£'((5yfuturc|a;T)/5(a;T|<5ypast)] ' 

(3.4) 

where p{xr\Sypast) is the analog of the filtering proba- 
bihty density P(xr |(5j/past) and (future |a;r) is the ana- 
log of the retrodictive likelihood function P^Sytuturclxr)- 
One can then solve for the density and effect operators 
separately, before combining them to form the classical 
smoothing probability density. 



C. Filtering 

Since the hybrid density operator can be regarded as 
a special case of the density operator, the same tools in 
quantum measurement theory can be used to derive a 
filtering equation for the hybrid operator. First, write 
p{xt+st\Syito,t]) in terms of p{xt\Sy[to,t]) as 



Because a quantum system is involved, one may be 
tempted to use a hybrid density operator d, 0, [3, [i3| 
to represent one's knowledge about the hybrid classical- 
quantum system. The hybrid density operator p{xr) de- 
scribes the joint classical and quantum statistics of a 
hybrid system, with the marginal classical probability 
density for Xr and the marginal density operator for the 
quantum system given by 

P{Xr)=tT[p{Xr)], (3.1) 
P{t) = / dXrp{Xr), (3.2) 



p{xt+5t\5y[toA) ^ j dxtJC[xt+st\xt)p{xt\5y[ta.t]), (3.5) 

where /C is a completely positive map that governs the 
Markovian evolution of the hybrid state independent of 
the measurement process. Equation p.5p may be re- 
garded as a quantum version of the classical Chapman- 
Kolmogorov equation. For infinitesimal 5t^ 

j dxtlC{xt+st \xt)p{xt) ~ [(i + StC) p{xt x)] ^^^^^^^ ■ 

(3.6) 
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The hybrid superoperator C can be expressed as 



Cp{x) = Cop{x) + Ci{x)p{x) 



d 



^ dxu 



2 ^ dxndxi, 

fi.y ^ 



(3.7) 



where Cq governs the evolution of the quantum system, 
Li governs the couphng of to the quantum system, via 
an interaction Hamiltonian for example, and the last two 
terms governs the classical evolution of Xt ■ 

Next, write p{xt\Syito^t]) in terms of p{xt\Sy[to,t-St]) us- 
ing the quantum Bayes theorem [18| as 



Kxt\6y[t„,t]) = pixt\Sy[t„^t-st],Syt) 

_ J{Syt\xt)p{xt\Sy[to^t-5t]) 
J dxt tr (numerator) 



(3.8) 



The measurement superoperator J'{Syt\xt), a quantum 
version of P{6yt\xt), is defined as 

Ji5yt\xt)p{xt) = M{6yt\xt)p{xt)M\Syt\xt). (3.9) 

For infinitesimal 6t and measurements with Gaussian 
noise, the measurement operator M can be approximated 

as \m 



M{Szt\xt) oc i + ^7^(t) 



-c^{xt,t)5z^t 



St 



—cl{xt,t)Cf,{xt,t) 



(3.10) 



where Szt is a vectoral observation process, c{xt,t) is a 
vector of hybrid operators, generalized from the purely 
quantum c operators in Ref. [Sj so that the observations 
may also depend directly on the classical degrees of free- 
dom, and 7^(<) is assumed to be positive. To cast the 
theory in a form similar to the classical one, perform uni- 
tary transformations on 5zt and c. 



5yt = U5zt, 
C{xt,t) = Uc{xt,t), 



(3.11) 
(3.12) 



where [/ is a unitary matrix, and rewrite the measure- 
ment operator as 

M{5yt\xt) ^i + \c^{xut)R~\t)5yt 



St 



C^^ixt,t)R-\t)C{xt,t). (3.13) 



C{xt, t) is a generalization of C(a;t, t) in the classical case, 
and R{t) is again a positive-definite matrix that charac- 
terizes the observation uncertainties and is real and sym- 
metric with eigenvalues 1/7^^. Note that 'I' is defined as 
the adjoint of each vector element, and is defined as 
the matrix transpose of the vector. For example. 



c't^i?-ic'^^(7t(i?-i)^,a 



(3.14) 



The evolution of p{xt\Sy[tg^t-St]) can thus be calculated 
by iterating the formula 

p{xt+st\Syito.t]) 

, ^. I ^J{Syt\xtyp{xt\8y\^ta,t-8t\) , 

dxtK.[xt+f,t\xt) TT- — — 7 — - — . (3.15) 

J axt tr(numerator) 

Taking the continuous time limit via Ito calculus and 
defining the conditional hybrid density operator at time 
t as 



F{x,t) = p{xt = x\dy[to,t)), 



(3.16) 



one obtains [1] 



dF = dtCF 
dt 

1 

+ 2 



2C^R-^FC^ - C^'^R'^^CF ~ FC'^'^R-^C 



C-{C)p) R-^dr^tF + Yi.c. 



where 



dx tr 



C{x,t)F{x,t) 



dt 



d^t = dyt--{C + C^) 



(3.17) 

(3.18) 
(3.19) 



is a Wiener increment with covariance matrix R{t)dt [l9[, 
H.c. denotes the Hermitian conjugate, and the initial 
condition is the a priori hybrid density operator p{xtQ ) ■ 
Equation ()3.17p is a quantum version of the KS equa- 
tion (|2.43p and can be regarded as a^ecial case of the 
Belavkin quantum filtering equation [20] . 

A linear version of the KS equation for an unnormal- 
ized F{x, t) is 



df = dtCf 
dt 
~8 



J (^C^R-^fC"^ - C'^'^R-^Cf ~ fC^^R-^c) 
+ ^[C' R-Uytf + ^.c}j , (3.20) 

and the normalization is 



F{x,t) 



J dxti[f{x,t)]' 



(3.21) 



Equation (|3.20p is a quantum generalization of the DMZ 
equation (|2.48|) . 



D. Retrodiction and smoothing 



Taking a similar approach to the one in Sec. Ill Dl and 
using the quantum regression theorem, one can express 
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the future observation statistics as [21| 



PiSytutmc) = / dxr tr 



dxT tr 
J{dyT~st\xT-5t 



E{Sy{utui-c\xT)p{xT) (3.22) 

J dxT-StlC{xT\xT-St) 



J dxT~25tK,[xT-St \xT-2St) 
J{6yT-2St\xT-2St) ■ ■ ■ 



dXrlC{Xr+St\Xr)J{SyT\xT)p{xT) , 

(3.23) 



which are analogous to Eq. (|2.50p and Eq. (|2.58p . re- 
spectively. Comparing Eq. (|3.22[) with Eq. (|3.23p . and 
defining the adjoint of a superoperator O as O*, such 
that 



tr 



E{x)Op{x) =tr| 0*E{x) /5(x)| , 
the hybrid effect operator can be written as 

EidytuturclXr) 
= J*{Syr\Xr) / dXr+StIC* {Xr+St\Xr) ■ ■ ■ 



(3.24) 



J*{SyT-25t\xT-2St) J dxT~StK.*{xT~St\xT-2St) 

J*{SyT-st\xT-st) / dxTlC*{xT\xT-st)i- (3.25) 



The operation K.* = J dx' K.* {x' \x) ■ may also be regarded 
as a hybrid superoperator on a hybrid operator, and is 
the adjoint of /C = J dx'IC{x\x')-, defined by 



E{x), lCp{x)J = {K.*E{x),p{x)'j , (3.26) 
with respect to the Hilbert- Schmidt inner product 



E{x),p{x))= / dxtr E{x)p{x) 



(3.27) 



One can then rewrite Eqs. ((3?22)) . ((3?23)) . and ((3?25)) more 
elegantly as 

(e{x),p{x)) = {l,]CJ... ICJp{x)) , (3.28) 
E{x) = J*IC* ...J*IC*i. (3.29) 

In the continuous time limit, a linear stochastic dif- 
ferential equation for the unnormalized effect operator 
q(x,t) oc E{dy[t T)\xt — x) can be derived. The result is 



—dg ~ dtC* g 

'^^ ^-'gC - gC^^R-'C - C^^R-'Cg 



+ ^ (2C^TR-'gC 
+ i (ffC'^i?-idyt+H.c.) , 



(3.30) 



to be solved backward in time in the backward Ito sense, 
with the final condition 



g{x,t) cx i. 



(3.31) 



Equation (|3.30p is the adjoint equation of the forward 
quantum DMZ equation p.20p with respect to the inner 
product defined by Eq. p.27p . It is a generalization of 
the classical backward DMZ equation (|2.6ip . 

Finally, after solving Eq. (|3.20p for f{x, r) and 
Eq. p.30p for g{x,T), the smoothing probability density 



h{x,T) = P{xr = x\dy[t„,T)) 



ti-[g(a^,T)/(x,r)] 

/dx tr[g(x,r)/(a;,T)] 

(3.32) 



The denominator of Eq. (I3.32p ensures that h{x,T) is 
normalized, so /(x, t) and g{x, t) need not be normalized 
separately. Table U lists some important quantities in 
classical smoothing with their generalizations in hybrid 
smoothing for comparison. 



E. Smoothing in terms of Wigner distributions 



To solve Eqs. ([3:201) . dSSQl), and dS^S]) . one way is to 
convert them to equations for quasiprobability distribu- 
tions [l^ . The Wigner distribution is especially useful for 
quantum systems with continuous degrees of freedom. It 
is defined as [22, 23] 



/(9,P) = ^ 



du (q 

2 



/ 



<? + o ) exp [ip 



(3.33) 



where q and p are normalized position and momentum 
vectors. It has the desirable property 



dqdpg{q,p),f{q,p) 



1 

2^ 



■ tr 



(3.34) 



which is unique among generalized quasiprobability dis- 
tributions 23]. The smoothing probability density given 
by Eq. p.32p can then be rewritten as 



h(x, t) 



J dqdpg{q,p,x,T)f{q,p,x,T) 
J dqdpdxg{q, p, x, T)f{q, p, x, r) ' 



(3.35) 



where f{q,p,x,T) and g{q,p,x,T) are the Wigner distri- 
butions of / and g, respectively. Equation (|3.35p resem- 
bles the classical expression (|2.68l) with the quantum de- 
grees of freedom q and p marginalized. If f{q,p,x,tQ) is 
nonnegative and the stochastic equations for f(q^p,x,t) 
and g{q,p,x,t) converted from Eqs. (|3.20p and (|3.30p 
have the same form as the classical DMZ equations given 
by Eqs. (|2.48p and (|2.6ip . the hybrid smoothing problem 
becomes equivalent to a classical one and can be solved 
using well known classical smoothers. For example, if 
f{q,p,x,t) and g{q,p,x,t) are Gaussian, h{x,T) is also 
Gaussian, and their means and covariances can be solved 
using the linear MFP smoother described in Sec. Ill El 
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Classical 


Description 


Hybrid 


Description 


Til ^ 1 ™ Si \ 

F(xt+st\xt,(iyt) 


transition probability density, ap- 
pears in Chapman-Kolmogorov equa- 
tion (I2.36P 


fi(Xt+St\Xt) 


transition superoperator, appears in 
quantum Chapman-Kolmogorov equa- 
tion 


P{Syt\xt) 


observation probability density, ap- 
pears in Baves theorem (|2.39p 


J{3yt\xt) 


measurement superoperator, appears 
in quantum Bayes theorem (13.81) 


P[Xt\dyito,t)), F{x,t) 


filtering probability density, obeys 
Kusliner-Stratonovich eciuation {12.43P 


P{xt\dy[to,t)), F{x,t) 


filtering hybrid density operator, 
obevs Belavkin eouation (|3.17p 


f{x,t) 


unnormalized F{x,t), obeys Duncan- 
Mortensen-Zakai (DMZ) equation 
(12.481) 


fix,t) 


unnormalized f{x,t), obeys quantum 
DMZ equation (|3.20|) 


P{dyyt,T)\xt) 


retrodictive likelihood function 


E{dyit,T)\xt) 


hybrid effect operator 


g{x,t) 


unnormalized P[dy^f^j')\xt)^ obeys 
backward DMZ equation (|2.6l|l 


g{x,t) 


unnormalized E[dy^f^x)\xt)^ obeys 
backward quantum DMZ equation 
(|3.30|) 


P{Xr\dyito,T)), h{x,T) 


smoothing probability density, obeys 
Eq. ((m)) 


P{xr\dyito,T)), h{x,T) 


smoothing probability density, obeys 
Eq. (I3.32P 



TABLE I: Important quantities in classical smoothing and their generalizations in hybrid smoothing. 



IV. PHASE-LOCKED LOOP DESIGN FOR 
NARROWBAND SQUEEZED BEAMS 



Consider the PLL setup depicted in Fig. [4l The opti- 
cal parametric oscillator (OPO) produces a squeezed vac- 
uum with a squeezed p quadrature and an antisqueezed 
q quadrature. The squeezed vacuum is then displaced by 
a real constant b to produce a phase-squeezed beam, the 
phase of which is modulated by (j)t = xu, an element of 
the vectoral random process xt described by the system 
Ito equation (j2.2p . The output beam is measured con- 
tinuously by a homodyne PLL, and the local-oscillator 
phase (jj't is continuously updated according to the real- 
time measurement record. 

The use of PLL for phase estimation in the presence of 
quantum noise has been mentioned as far back as 1971 
by Personick [24| . Wiseman suggested an adaptive homo- 
dyne scheme to measure a constant phase [25j , which was 
then experimentally demonstrated by Armen et al. for 
the optical coherent state [1^. Berry and Wiseman [27| 
and Pope et al. [11] studied the problem with </>( being 
a Wiener process. Berry and Wiseman later generalized 
the theory to account for narrowband squeezed beams 
Q. Tsang et al. also studied the problem for the case of 
Xt being a Gaussian proces s [29l.l30j . but the squeezing 



model considered in Refs. 23? 30jis not realistic. Us- 
ing the hybrid smoothing theory developed in Sec. IIIIl 
one can now generalize these earlier results to the case 
of an arbitrary diffusive Markov process and a realistic 
squeezing model. 

Let p{xt) be the hybrid density operator for the com- 
bined quantum-OPO-classical-modulator system. The 
evolution of the OPO below threshold in the interaction 



picture is governed by 
CqP{x) = 



Ha,p{x) 



^ {qp + pq) , 



(4.1) 
(4.2) 
(4.3) 



where a is the annihilation operator for the cavity optical 
mode, and q and p are the antisqueezed and squeezed 
quadrature operators, respectively, defined as 

q = — (4.4) 



V2 



P ^ 



V2i 

with the commutation relation 
[q,p] = i 



(4.5) 



(4.6) 



The classical phase modulator does not influence the evo- 
lution of the OPO, so 

Ci - 0, (4.7) 

but it modulates the OPO output. C{xt,t) in this case 
is 

Cixt,t) = -2i {b + y^a) exp (if/jt - i<\>'t) , (4.8) 

where 7 is the transmission coefficient of the partially 
reflecting OPO output mirror, i? = 1, and the sym- 
bol and sign conventions here roughly follows those of 
Refs. [2^, [so]. To ensure the correct unconditional quan- 
tum dynamics, the Hamiltonian should be changed to 
(Ref. [ll. Sec. 11.4.3) 

Hl, = Ho- i'-^^ia - at), CoK^) = [^^, /5(x) 

(4.9) 
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FIG. 4: (Color online). Homodyne phase-locked loop (PLL) for phase estimation with a narrowband squeezed optical beam 
produced from an optical parametric oscillator (OPO). 



in order to eliminate the spurious effect of the displace- 
ment term in C on the OPO. After some algebra, the 
forward stochastic equation for the Wigner distribution 
f{q,p,x,t) becomes 



d 



1 ^ 92 



2 dxndxy 
-X - J) I iPf) 



) f 



1 

4 

dyt 
cos( 



5V 
sin(( 



7 d 
2 dp 



/• 



(4.10) 



This is precisely the classical DMZ equation (12.48^ with 
correlated system and observation noises. The equivalent 
classical system equations are then 



dqt= ( X - 2 ) ltdt 



dpt = (-X - I) Ptdt + ^ '^dPt, 
dxt = A{xt,t)dt + B{xt,t)dWt, 
and the equivalent observation equation is 



(4.11) 



dyt = 26sin( 
dCt 



sm 



,-^'t)dt + dCt, 

4''t) {y^qtdt - dat 
-0;) [^ptdt-d(3t 



(4.12) 



where dat and d(3t are independent Wiener increments 
with covariance dt. dat and d/^t, which appear in both 
the system equation and the observation equation, are 
simply quadratures of the vacuum field, coupled to both 
the cavity mode and the output field via the OPO output 
mirror. Equations (|4.1ip and (|4.12p coincide with the 
model of Berry and Wiseman in Ref. Q when xt is a 



Wiener process, and Eq. (|4.10p is the continuous limit of 
their approach to phase estimation. This approach can 
also be regarded as an example of the general method of 
accounting for colored observation noise by modeling the 
noise as part of the system 0, H, 0] . 

If X = 0, dC,t/dt is an additive white Gaussian noise, 
and the model is reduced to that studied in Refs. J^, 
[29I, [S^l . In that case, it is desirable to make 4>'t follow (j)t 

can be approximated 



as closely as possible, so that dyt 
as 

dyt~2b{(j)t-(i)'t)dt + dCt, 



(4.13) 



and the Kalman filter can be used if xt is Gaussian [30|. 
Provided that Eq. (|4.13p is valid, one should make 4>'t the 
conditional expectation of (pt = xu, given by 



(4.14) 



0t = (0t)/ = / dqdpdxxif{q,p,x,t). 



For phase-squeezed beams, it also seems desirable to 
make 0j close to (j)t in order to minimize the magnitude 
of dCf Equation (|4T4)) may not provide the optimal 0^ 
in general, however, as it does not necessarily minimize 
the magnitude of d(^t or the estimation errors. The opti- 
mal control law for 0j should be studied in the context 
of control theory. 

While (f)f needs to be updated in real time and must 
be calculated via filtering, the estimation accuracy can 
be improved by smoothing. The backward DMZ equa- 
tion for g{q,p,x,t) is the adjoint equation with respect 
to Eq. (|4.10p , given by 



-d9 = dt{j:A,^ + lT.iBQBn,. 



dxadxv 



1\ dg_ 
dp 



1 fd^g 9^9 



+ dyt 



sm 



cos 



2b 

2-ip- 



279- 

7 d 
28^ 



7 d 



(4.15) 
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and the smoothing probabihty density h{x, r) is given by 
Eq. (|3.35p . The use of hnear smoothing for the case of Xt 
being a Gaussian process and d(t / dt being a white Gaus- 
sian noise has been studied in Refs. [1^, [sOj- Practical 
strategies of solving Eqs. (|4.10p and (|4.15p in general are 
beyond the scope of this paper, but classical nonlinear fil- 
tering and smoothing techniques should help P, [3, 0, 0] ■ 
One can also use the hybrid smoothing theory to study 
the general problem of force estimation via a squeezed 
probe beam and a homodyne PLL, by modeling the phase 
modulator as a quantum mechanical oscillator instead 
and combining the problem studied in this section with 
the force estimation problem studied in Ref. [5|]. 



WEAK VALUES AS QUANTUM 
SMOOTHING ESTIMATES 



Smoothing Estimation 



Measurements 

M{dyt) ■ 



Quantum System 


dyt 


Observations 





FIG. 5: (Color online). The quantum smoothing problem. 

Previous sections focus on the estimation of classical 
signals, but there is no reason why one cannot apply 
smoothing to quantum degrees of freedom as well, as 
shown in Fig. \5\ First consider the predicted density 
operator at time r conditioned upon past observations, 
given by 



Pir) 



fir) 



(5.1) 



where the classical degrees of freedom are neglected for 
simplicity. The predicted expectation of an observable, 
such as the position of a quantum mechanical oscillator, 
is 



(6);^tr 6 Pir) 



tr[0./(r)] 
tr[/(r)] ' 



(5.2) 



One may also use retrodiction, after some measurements 
of a quantum system have been made, to estimate its 
initial quantum state before the measurements (sil . [s^ , 
using the retrodictive density operator defined as 



Aot(T) = 



9ir) 



trL9(r)J 

The retrodicted expectation of an observable is 

tr[g(r)0] 



;(0) =tr /5,et (t)0 



tr . 



0] 



(5.3) 



(5.4) 



Causality prevents one from going back in time to verify 
the retrodicted expectation, but if the degree of freedom 
with respect to O at time r is entangled with another 
"probe" system, then one can verify the retrodicted ex- 
pectation by measuring the probe and inferring O [33 | . 

The idea of verifying retrodiction by entangling the 
system at time r with a probe can also be extended to 
the case of smoothing, as proposed by Aharonov et al. 
Q. In the middle of a sequence of measurements, if one 
weakly couples the system to a probe for a short time, so 
that the system is weakly entangled with the probe, and 
the probe is subsequently measured, the measurement 
outcome on average can be characterized by the so-called 
weak value of an observable, defined as [1, [s^l 



f 



tr[g(r)0/(T)] 
tr[g(r)/(r)] ' 



(5.5) 



The weak value becomes a prediction given by Eq. (|5.2p 
when future observations are neglected, such that ^(t) = 
1, and becomes a retrodiction given by Eq. (j5.4[) when 
past observations are neglected and there is no a priori 
information about the quantum system at time t, such 
that /(t) = 1. When /(r) and ^(r) are incoherent mix- 
tures of O eigenstates. 



fir)^Y.fiO,r)\0)(0\, 
o 

SM-E5(0,r)|0)(0|, 



the weak value becomes 



f 



Eo09iO:r)f{0,T) 



(5.6) 
(5.7) 

(5.8) 



and is consistent with the classical time-symmetric 
smoothing theory described in Sec. [TTl Hence, the weak 
value can be regarded as a quantum generalization of the 
smoothing estimate, conditioned upon past and future 
observations. 

One can also establish a correspondence between a 
classical theory and a quantum theory via quasiprobabil- 
ity distributions. Given the smoothing probability den- 
sity in terms of the Wigner distributions in Eq. (j3.35p . 
one may be tempted to undo the marginalizations over 
the quantum degrees of freedom and define a smoothing 
quasiprobability distribution as 



h{q,p,T) 



9i<l,P, T)f{q,p,T) 
J dqdpg{q,p, T)f{q,p, 



(5.9) 



where f{q,p,T) and g{q,p,T) are the Wigner distri- 
butions of /(r) and ^(t), respectively. Intriguingly, 
h[q,p,T), being the product of two Wigner distributions, 
can exhibit quantum position and momentum uncertain- 
ties that violate the Heisenberg uncertainty principle. 
This has been shown in Ref. j30j , when the position of a 
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quantum mechanical oscillator is monitored via contin- 
uous measurements and smoothing is applied to the ob- 
servations. From the perspective of classical estimation 
theory, it is perhaps not surprising that smoothing can 
improve upon an uncertainty relation based on a predic- 
tive theory. The important question is whether the sub- 
Heisenberg uncertainties can be verified experimentally. 
Ref. [30] argues that it can be done only by Bayesian 
estimation, but in the following I shall propose another 
method based on weak measurements. 

It can be shown that the expectation of q using 
h{q,p,T) is 



{q)h= I dqdpqh{q,p,T) 
tr[.9(r)g/(T)] 



Re- 



tr[5(r)/(r)] 



R.eg(g);, 



(5.10) 
(5.11) 



which is the real part of the weak value, and likewise for 
so the smoothing position and momentum estimates 
are closely related to their weak values. More generally, 
consider the joint probability density for a quantum po- 
sition measurement followed by a quantum momentum 
measurement, conditioned upon past and future obser- 
vations: 



From the perspective of classical probability theory, 
Eq. (|5.16p can be interpreted as the probability density of 
noisy position and momentum measurements with noise 
variances l/cq and 1/ep, when the measured object has 
a classical phase-space density given by P{q,p). In the 
limit of infinitesimally weak measurements, eq,ep 0, 
and 



lim P{q,p) 



h{q,p,T) 



(5.18) 



Thus, h(q^p,T) can be obtained approximately from an 
experiment with small eg and ep by measuring P{yq,yp) 
for the same g and / and deconvolving Eq. (|5.16p . In 
practice, Cg and Cp only need to be small enough such 
that P{q,p) ~ h{q,p, r). This allows one, at least in prin- 
ciple, to experimentally demonstrate the sub-Heisenberg 
uncertainties predicted in Ref. [13] in a frequentist way, 
not just by Bayesian estimation as described in Ref. |3(j ]. 
Note, however, that h{q,p,T) can still go negative, so 
it cannot always be regarded as a classical probability 
density. This underlines the wave nature of a quantum 
object and may be related to the negative probabilities 
encountered in the use of weak values to explain Hardy's 
paradox [s^ . 



1 



Piy„yp) = ^tr [g{T)A%{yp)Mg{yg)f{T)Ml{yg)M;{yp) 

(5.12) 



(5.13) 



C EE y dygdyptr \g{T)Mp{yp)Mg{yg 
X j{T)Ml{yg)Ml{yp 
where the measurement operators 



exp 



^{yg-qf \q){ql 



(5.14) 



Mp{yp) = / dp (1^) " exp - j(2/p 



p? b>(p| 



(5.15) 



are assumed to be Gaussian and backaction evading. Af- 
ter some algebra. 



pto„»,)^/<i,*(U)'(a- 



X exp 



(5.16) 



P{q.p) 



— — / dudv exp 



u 
2 



9{t) P - I) cxp(iwg) 



X {q 



fir) 



q+ ■^)exp(ipu). 



(5.17) 



VI. CONCLUSION 

In conclusion, I have used a discrete-time approach 
to derive the classical and quantum theories of time- 
symmetric smoothing. The hybrid smoothing theory 
is applied to the design of PLL, and the relation be- 
tween the proposed theory and Aharonov et aL's weak 
value theory is discussed. Possible generalizations of the 
theory include taking jumps into account for the clas- 
sical random process t9| and adding quantum measure- 
ments with Poisson statistics, such as photon counting 
[Tsl . [2ll . m, H^. Potential applications not discussed 
in this paper include cavity quantum electrodynamics 
fT8l[2ll.[2^.l23t. photodetection theory 16, 18, 23], atomic 
magnetometry [35| , and quantum information processing 
in general. On a more fundamental level, it might also be 
interesting to generalize the weak value theory and the 
smoothing quasiprobability distribution to other kinds of 
quantum degrees of freedom in addition to position and 
momentum, such as spin, photon number, and phase. A 
general quantum smoothing theory would complete the 
correspondence between classical and quantum estima- 
tion theories. 
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